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ABSTRACT 

This paper presents the results of searches for anisotropy in arrival directions of ultra-high energy 
cosmic rays detected with the Yakutsk Array during the 1974-2008 observational period together with 
available data from other giant extensive air shower arrays working at present. A method of analysis 
based on a comparison of the minimal width of distributions in equatorial coordinates is applied. As 
a result, a hypothesis of isotropy in arrival directions is rejected at the 99.5% significance level. The 
observed decrease in the minimal width of distribution can be explained by the presence of cosmic 
ray sources in energy intervals and sky regions according to the recent indications inferred from data 
of the Yakutsk Array and Telescope Array experiments. 

Subject headings: astroparticle physics - cosmic rays 


1. INTRODUCTION 

The origin of ultra-high energy cosmic rays (UHECRs) 
is a long-standing challenge in cosmic ray (CR) physics. 
Extensive air shower (EAS) arrays detecting CRs at en¬ 
ergies above 1 EeV (=10^® eV) observe mainly isotropic 
arrival directions with no sign of fluxes from s ources sig¬ 
nificant l y exceeding instrumental erro rs (e.s.. lAab et al.l 
(|2014all : iKampert fc TinvakovI (I2014IB . 

At the same time, there are indications of the small 
size anisotropy in arrival directions revealed by means 
of a comp arison of CR intens it ies in adjacent sky re- 
:ions f e.g.. Ilvanov et al.l (l2003ll : ISommers fc Westerhofll 
20091) : lAbbasi et arr ( 2014l B. To confirm or refute the in¬ 
dications it is useful to diversify approaches and methods 
in analysis of CR arrival directions besides independent 
experimental data. 

Recently, a method of the minimal width of arrival 
direction distribution (MWADD) was used^ to search for 
the sources of UHE CRs in a circle of right ascensions 
(jivanov et ahll^lhaf) . The method is a specific variant 
of testing for the equality of variances of p opulations, a s 
an alternative to the equality of the means (IIvanovll2013D . 

In this paper, a two-dimensional generalized method 
is developed for application in equatorial coordinates to 
test a null hypothesis, Hq, of isotropy in the arrival di¬ 
rections of CRs and an alternative hypothesis with fitted 
position and luminosity of a possible source of CRs. Our 
aim is to reject or confirm and refine characteristics of 
the possible CR sources indicated in the previous papers. 

The present enhancement of the method consists in 
consideration of the non-uniformity of the array accep¬ 
tance area in CR arrival angles due to the unequal time- 
integrated flux from different directions of the sky. 


2. THE YAKUTSK ARRAY EXPERIMENT AND 
SAMPLING OF THE DATA SET 

The main purpose of the Yakutsk Array^ is to in¬ 
vestigate CRs measuring EAS in the energy range of 


^ for a method in the context of directional statistics see Ap¬ 
pendix 

^ Website: http://eas.ysn.ru 


10^® — I0^° eV. Construction of the array near Yakutsk, 
Russia, at geographical coordinates 61.7°N, 129.4°i/, 105 
m above sea level (102 0 g/cm^), was completed in 1973 
(jDvakonov et al.l|l991h . During years, the array has been 
reconfigured several times. Before 1990, the total area 
covered by detectors was ^17 km^; now, it is 8.2 km^. 
At present, it consists of 58 ground-based and four un¬ 
derground scintillation counters to measure charged par¬ 
ticles (electrons and muons), and 48 detectors of the 
air C herenkov light (|Egorova et al.l l2001t Ilvanov et al.l 

[2nTnl) . 

EAS events selection from the background is realized 
with a two-level trigger of detector signals: The first level 
is a coincidence of signals from two scintillation counters 
in a station within 2 /is; the second level is a coincidence 
of signals from at least three nearby stations within 40 fis. 
Functioning pro cedures and the types o f array d e tector s 
are described in iDvakonov et al.l (jl991h : llvanovi (|2009D : 
Ilvanov et al.l (120131) . 

The location of the shower core is based on the fitting 
of the particle lateral distribution by the Greisen-type 
trial function. Core location errors are ^ 30 to ^ 50 m 
depending mai nly on the number of tr iggered stations in 
the EAS event (|I)vakonov et al.iri99lD . 

Arrival angles of EAS primary particles are calculated 
in the plane shower front approximation using the trigger 
times of stations. A clock pulse transmitter at the center 
of the array provides a pulse timing of ^ 100 ns accuracy. 
Errors in arrival angles depend on the primary energy, 
decreasing from ^ 7° at U = 1 EeV to ^ 3 ° above E = 
10 EeV. For more detailed information see Ilvanov et al.l 
(l2009l[2(rm[) . 

In contrast to two other giant arrays function- 
ing at presen t, Pie rre Auger Observatory (PAO, 
Abraham et al.l (l2004l) ) and the Telescope Array (TA, 
Kawai et al.l (j2009t )) with fluorescent detectors, atmo¬ 
spheric Cherenkov light is used here to esti mate the en¬ 
ergy of primary p articles initiating EAS (jivanov et al.l 

[2(M[2(m[2nT5H) . 

In this work, the same sample of the Yakutsk Ar- 
ray data is used as in the previous paper (jivanov et al.l 
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Fig. 1.— Expected declination distribution of the isotropic cos¬ 
mic rays detected with EAS arrays. N = 10®. 


Fig. 2.— Half width of the isotropic arrival direction distribution 
as a function of declination. N = 10®. 


consisting of EAS events detected in the period 
January 1974 - June 2008, with shower axes within the 
array area with zenith angles 0 < 60°. The unified energy 
estimation procedure is applied to all showers through¬ 
out the period, following lEgorova et al.l (l2001h . An ad¬ 
ditional cut is caused by the increased threshold energy, 
which will be set out in the next subsection. 


2.1. Exposure of the ground EAS array for celestial 

regions 

Because the array aperture is bounded in the zenith 
angle interval 9 G {^,0thr), CRs can be detected only 
from a part of the sky in a horizontal system. A diurnal 
cycle of the array functioning provides nearly uniform ex¬ 
posure in right ascensions and apparent non-uniformity 
in declinations. 

A simple and convenient way to calculate the non- 
uniform exposure of the array is to use a Monte Carlo 
MC) n iethod (for basics see, e.g., iMetropolis fc UlamI 
194911 : iHammerslev &: HandscomU (|1975ll : an appli¬ 
cation to t h e array ex posure is demonstrated by 
llvanov et al.l (|1997l 1200,11 11. An algorithm is based on 
the angular distribution of isotropic rays in a horizon¬ 
tal system attached to the flat array on the ground. 
Within the infinitesimal time interval {t,t dt), with 
stationary Earth, the azimuthal distribution is uniform 
and the zenith angle distribution is formed as follows: 
if {9 < 9thr) then fi{0) = sin(26')/(l - cos'^ 9thr), else 
fi{9) = 0. If we assume the inverse time t —>■ —t, so that 
all the exposed rays move from the array to the sky, then 
the rays are integrated over the diurnal cycle due to the 
Earth reverse rotation. 

In this algorithm, random directions, {(j)i,9i),i = 
I,..., N, are sampled in the horizontal system from a uni¬ 
form distribution in the azimuth and fi{9) in zenith an¬ 
gles. A uniform distribution of the sidereal time is used 
as well. Then directions are converted, f or instance, t o 
equator i al an gles {(t>i,9i) —>• {ai,5i) fe.e.. iGreeiil (fM) : 
ICollinsI (j2004ll l to form the expected-for-isotropy distri¬ 
bution of CR arrival directions on the celestial sphere. 

Resultant declination distributions for the ground ar¬ 
rays working at present are shown in Fig. [U while the 
right ascension distribution is uniform. Actually, the ob¬ 
tained distributions are formed by the directional expo¬ 
sures of EAS arrays, that is, the effective time-integrated 
collecting area for a flux from each direction of the sky. 

Minor deviations from the uniform distribution in right 
ascension are caused by the diurnal and seasonal varia¬ 


tions of the array exposure. The efficiency of array de¬ 
tection is affected by the weather effects, shutdowns of 
detectors, the ge omagnetic modulation of the EAS event 
rate, and so on (iPravdin et al.l[20MI : llvanov et al.lll999l 
|2 oo 3). Attenuation of showers in the atmosphere results 
in a deviation from the ‘isotropic’ zenith angle distribu¬ 
tion, /i, at low energies. However, above some threshold 
energy, Et^r, all these effects can be neglected against 
statistical errors rising with energy. In this paper, the fol¬ 
lowing values are c hosen: Efhr = 3.2 Ee V, 9thr = 60° for 
the Yakutsk Array (llvanov et al.l 2015a l. Ethr = 52 EeV, 
9thr = 80° for PAO (lAab et~aiTl2014bli. and Efh r = 57 
EeV, 0thr = 55° for TA data ( Abbasi et al.ll2014 1. 


3. ANALYSIS OF CR ARRIVAL DIRECTIONS IN 
EQUATORIAL COORDINATES 

3.1. A method of the minimal width of distribution to 
analyze arrival directions of CRs 

Isotropic distribution of arrival directions has no mean 
value because the limits of the region on the sphere that 
should be integrated over are undefined. Meanwhile, if 
we assume an arbitrary direction, (a, J), as a trial mean, 
then we can hnd a dispersion of the isotropic distribution, 
namely, the width 2a;i, which is independent of the trial 
mean: 

^ p2t7 pT7 ^ 

zji = — d(j) ifsinipdip = -, 

47r Jo Jo 

where if is the angular distance to (a, 6). In the case of 
N data points on the celestial sphere, a sum of angular 
distances is applicable instead: 

1 ^ 
i=0 

where the asymptotic limit is equal to tt/2 as N ap¬ 
proaches oo. 

So, the width of the isotropic distribution on the sphere 
is 2u!i = 180°. On the other hand, if there is a source of 
CRs with the angular size S' <C Wi lurking in an isotropic 
background, then the aggregate width of the distribution: 
i) reaches the minimum when the trial mean points to the 
source; ii) has a minimum which is distinctly less than 
2u}i, depending on the fraction of the source luminosity 
in the overall flux of CRs. 

For instance, if the flux from the source is half of the 
total, then ujmin = 45°, while for a fraction of 0.1, half of 
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the minimal width is uJmin = 81°. It seems that by mea¬ 
suring the width of the distribution of arrival directions, 
one is able to reject the null hypothesis and to find the 
coordinates and the fraction of CR flux from a source, if 
there is any. 

3.2. Application of the method to Monte Carlo data 

In this section, simulation results of the MWADD 
method applied to N points on the equatorial sphere 
taking into account the array exposure are given. In 
this case, the width of distribution is strongly influenced 
by the exposure, so the directional dependence uji{a,6) 
should be calculated for the particular array. 

It is straightforward to use the MC algorithm described 
above (Section 12.11) to compute the width with the trial 
mean scanning the whole a G (0, 360°),d G (—90°, 90°) 
equatorial area. The results for three arrays are shown 
in Fig. [5] At A > Ethr the distribution width is the 
same in right ascensions, so the width variation is shown 
for the trial mean scanning declinations. 

The method is applicable only in searching for a single 
source, SS, of CRs. Indeed, in the case of two sources 
located at the angular distance L from each other, with 
the fractions of CR fluxes / and (I — /), MWADD is 
ijj 2 = 2/(1 — f)L. For opposite sources with equal frac¬ 
tions, a half-width is ui 2 = 7’‘/2, just as in the isotropic 
alternative. In what follows, we will explicitly suppose 
an SS of CRs within a particular energy interval. 

The statistical power of the MWADD method is the 
efficiency depending on the sample size N. We have to 
find a lower limit of N needed to reject Hq at a confidence 
level of 99% when an alternative hypothesis. Hi, is true. 
To estimate Nmin, we used Hi consisting of SS as a d- 
function located in {ass,^ss), within the field of view of 
the Yakutsk Array, with the fraction of the total CR flux 
/, and an isotropic background which provides (1 — /) 
of the flux. The distribution width for each trial mean 
is normalized using the ‘exposed’ isotropic distribution 
width as a measure. 

The result of MC simulation is given in Fig. |3]in com¬ 
parison with the power of harmonic analysis in the right 
ascension^. The first harmonic amplitude, Ai, under Hi 
is a w eighed vector sum of 2 and 2I\/N (jivanov et al.l 
[inil^)- Using a probability P{> Ai) = exp{—NA\/A) = 
0.01, one can find Nmin for Hq as a function of /. A con¬ 
clusion to be drawn is that the minimal width method is 
more powerful than harmonic analysis in R.A. 

3.3. Application of the method to experimental data 

3.3.1. Testing the null hypothesis with the Yakutsk Array 

data 

The Yakutsk Array data at energies above 3.2 FeV 
are divided into four intervals with the widths A Ig if = 
0.25, 0.5. Below, energy scaling factors for FAS arrays 
deri ved from comparison of the observed energy s pec- 
tra (Ilvanovll2010l : ll)awson et al1l2013tlD.Ivanovll2014ll are 
used (PAO: 1.04; TA: 0.96; Yakutsk: 0.561). Scaled en¬ 
ergy is marked Ewg- 

The width of the observed distribution of arrival direc¬ 
tions is normalized using the expected isotropic distribu¬ 
tion width for a given trial mean. Only in the energy 

® In other words, the Rayleigh te st where th e first harmonic 
amplitude is a measure of dispersion Oupdl200in 
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Fig. 3.— Minimal sample size, Nmin, needed to reject isotropic 
distribution if there is a single source, SS, giving a fraction of the 
total CR flux. 


TABLE 1 

MC SIMULATION RESULTS: THE PROBABILITY, Pq, OE MWADD 
UNDER Ho TO BE LESS THAN OR EQUAL TO THE OBSERVED VALUE. 

The number of EAS events observed in energy bins, N^bs, 

AND A sample SIZE, M, USED IN SIMULATION ARE GIVEN EOR 
ARRAYS. 


Experiment, energy bin, EeV 

tYobs 

M 

Po,% 

PAO, Ewg > 54 

231 

10000 

57.3 

TA, > 55 

72 

100000 

0.1 

Yakutsk, (3.1, 5.6) 

939 

10000 

98.5 

Yakutsk, (5.6,10.0) 

285 

10000 

0.1 

Yakutsk, (10.0,17.7) 

95 

100000 

14.6 

Yakutsk, (17.7,56.1) 

42 

100000 

23.5 


interval Ewg G (5.6,10) FeV there is a definite mini¬ 
mum, iOobsl^i = 0.89, of the distribution width of CRs 
detected with the Yakutsk array, shown in Fig. |4l the 
central map. Another minimum of the width of arrival 
directions at energies above 55 Fe V is revealed i n dat a 
provided by the TA Collaboration (|Abbasi et al.l (l2014ll , 
mapped on the right). 

To estimate the probability of MWADD under H^ be¬ 
ing less than or equal to the observed value, the MC algo¬ 
rithm is used with the number of isotropic events in a set 
equal to the number of observed EAS events. Nobs, in the 
particular energy bin. The exact algorithm of MWADD 
calculation that was used for the data is then performed 
on the MC event set. 

The procedure is repeated M times to find the fraction 
of MC event sets where the minimal distribution width is 
equal or less than the experimental value. This fraction 
is interpreted as a probability to quote the significance 
of the anisotropy signal. The number of MC event sets 
used in simulation and the number of events in energy 
bins detected in experiments are presented in Table 1. 

The resultant probability for the Yakutsk array data in 
the energy bin Ewg G (5-6,10) FeV is Pq = 1.15 x 10“°, 
which is equivalent to ^ 3.1(7 deviation in the normal 
distribution terms. However, a penalty factor should be 
applied to the probability, which is calculated a posteri¬ 
ori. Assuming equally possible anisotropy in any of the 
four energy bins, with comparable deviations, one has a 
final probability P = 4.6 x 10“°, equivalent to ^ 2.6 ct. 
Consequently, the null hypothesis can be rejected basing 
on the Yakutsk Array data at the significance level of 
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99.5%. 

Observed valu es of MWADD ca lculated u sing available 
data from PAO (|Aab et al.ll2014bD and TA (|Abbasi et al.l 
I 2 OI 4 II are shown in comparison with the Yakutsk Array 
data in Fig. [5l There is no deviation from isotropic 
expectation in the data from PAO, while TA data exhibit 
a pronounced deviation of MWA DD in the energy r ange 
where a ‘hotspot’ was indicated (lAbbasi et al.ll20l3l . 

The probability of 72 isotropic EAS events above 
Ewg > 55 EeV having a MWADD less than that ob¬ 
served by the TA is Pq = 1-3 x 10“^ Str). A penalty 
factor can be calculated by the TA Collaboration only 
using all the data observed. 

3.3.2. Searching for the coordinates of possible CR sources 
and fraction of the flux 

Our alternative hypothesis, Hi, has free parameters to 
adjust to the observed MWADD: the position of SS in 
an equatorial system and the fraction of the total CR 
flux that has arrived from the source. The same itera¬ 
tive procedure is used to calculate the ‘most probable’ 
parameters yielding the observed distribution width. As 
an implementation, the MC program mentioned above is 
adapted to calculate MWADD under Hi with input free 
parameters. 

Fitted parameters are then used to calculate the ran¬ 
dom dispersion of MWADD for a fixed N equal to the 
number of detected CRs in a particular energy bin. By 
varying a parameter, its confidence interval is determined 
where the resultant deviation of MWADD is within ran¬ 
dom dispersion limits. 

Hypothesis Hi is applied to the Yakutsk Array data in 
the energy interval Ewg G (5.6,10) EeV and to TA data 
at Ewg > 55 EeV. The most probable parameters for 
the Yakutsk Array data are ay = 36° ^gg, fty = 48° ^^ 7 , 
fy = 0.11±0.08. The results for TA are apA = 144° ^ 30 ; 
^t^ = 42° ^23^ /r^ = 0 . 2 ± 0 . 1 . 

A hint of the possible source of CRs in the interval 
Ewg € (5.6,10) EeV, derived from the Yakutsk Array 
data u sing the MWADD m ethod in the right ascension 
circle (|Ivanov et al.ll201^ . is confirmed; the resulting 
ay intervals are within experimental errors. 

The coordinates of a hypothesized TA source are in 
agreement wit h that of a hotspo t revealed by the TA 
Collaboration (lAbbasi et al.ll20lll . An additional bonus 


in our case is an estimation of the probable fraction of 
CRs attributed to a source. 

4. CONCLUSION 

The Yakutsk Array data on arrival directions of CRs 
above 3.2 EeV in equatorial coordinates are analyzed us¬ 
ing the minimal width of distribution method. A pre¬ 
vious hint of large-scale anisotropy in the energy range 
5.6 < Ewg < 10 EeV is confirmed by the enhanced 
method. The null hypothesis is rejected at the 99.5% 
significance level. 

For comparison, our method of analysis is applied to 
available data from other giant EAS arrays. PAO data 
demonstrated no deviation of the MWADD from the 
isotropic distribution width at energies above 54 EeV. 
On the contrary, arrival directions of UHECRs detected 
with TA have a decreased minimal width at Ewg >55 
EeV with a statistical significance of ~ 3cr. 

The MWADD method is applied with an alternative 
hypothesis of a single CR source in a weighed combina¬ 
tion with the uniform background. Free parameters of 
a source are fitted using experimental data in the en¬ 
ergy intervals where Hq is ruled out: ay = 36° Igg, 
5y = 48° ^ 77 , fy = 0.11 ± 0.08 for the Yakutsk Array 
data in the energy interval Ewg S (5.6, 10) EeV, and 
apA = 144° 1^9^ Sta = 42° = o.2 ± 0.1 for TA 

data at Ewg > 55 EeV. 

Although our alternative hypothesis is not a unique ex¬ 
planation of the observed decrease in MWADD, the indi¬ 
cated coordinates and fitted CR fraction of the possible 
sources may be useful in a future search with enhanced 
statistics for anisotropy in arrival directions of UHECRs. 

The equatorial coordinates of the possible source de¬ 
rived from TA data are in agreement with a hotspot po¬ 
sition found by comparison of CR events summed within 
sky regions (lAbbasi et al.ll20l3l . Our approach is dif¬ 
ferent, being based on the overall distribution width of 
arrival directions rather than on the excess flux of CRs 
in a particular angular region. 

The author is grateful to the Yakutsk Array staff for 
the data acquisition and analysis. The work is supported 
in part by the Russian Academy of Sciences (Program 
10.2) and RFBR (grants 11-02-00158 and 13-02-12036). 


APPENDIX 


Our objectives in the paper are a circular uniform distribution on the unit sphere 5^ in ]R° and, as antithesis, a 
point source of CRs. For simplicity, all considerations in the Appendix will be illustrated on a circle, where i/'o is the 
point source position. 

The most common way in directional statistics is to use n-th moments of a distribution defined as = 
exp(m^i), where iV is a number of data points at ipi, with the mean angle ip = Argimi) and the circu¬ 
lar variance S' = 1 — |mi|. However, there are other measures of the distribution moments in use. For instance, we 
are using in the paper a measure of dispersion of angles X]i=i('^ ~ ~ l'0i ~ f/’oll) considered bv iMardle] (119721 1. In 

this approach the mean angle is that where dispersion reaches the minimum (for distributions under consideration the 
median coincides with the mean). 

The MWADD method is used in the paper to locate the mean and to test for uniformity of directions against 
alternative hypothesis with a point source. In Section 3.2 , the method is compared to the Rayleigh test known as one 
of the most powerful tests for unimodal data (|ju iaiMIl). 
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Fig. 5. — Minimal width of the distribution of arrival directions 
as a function of CR energy. A normalized difference in the widths 
of observed and isotropic distributions is derived using the data 
from EAS arrays. 












